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We study the mathematical character of the angular moment equations 
Oh' of radiative transfer in spherical symmetry and conclude that the system is 

6 

5h ' hyperbolic for general forms of the closure relation found in the literature. 

Hyperbolicity and causality preservation lead to mathematical conditions al- 
lowing to establish a useful characterization of the closure relations. We apply 
C3 ' numerical methods specifically designed to solve hyperbolic systems of con- 

servation laws (the so-called Godunov-type methods), to calculate numerical 
solutions of the radiation transport equations in a static background. The fea- 
sibility of the method in any kind of regime, from diffusion to free-streaming, 
is demonstrated by a number of numerical tests and the effect of the choice 
of the closure relation on the results is discussed. 
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2 J. A. Pons, J. Ibdnez and J. A. Miralles 
1 INTRODUCTION 

In standard problems of radiation^ hydrodynamics (RH) where radiation contains a large 
fraction of the energy and momentum density, the Boltzmann Equation (BE) must be cou- 
pled to the hydrodynamic equations in order to obtain the evolution of the system as well 
as the spectrum and angular distribution of the radiation field. However, an algorithm built 
to solve the BE numerically (a Boltzmann solver) in a non-stationary case is too time- 
consuming, from a computational point of view, to allow for the extension to more than one 
dimension of the existing numerical codes ( [Mezzacappa fc Bruenn 1993| ; [Yamada, Janka & 



Suzuki 1999|) . In many cases, instead of the BE, its angular moments are considered, obtain- 



ing then the multigroup (or multi-frequency) equations or the even simpler energy averaged 
equations. 

Standard RH methods flMihalas fc Mihalas 1984| ) have found their way into the litera- 



ture, and one can find several RH codes with different approaches for the radiation transport 



part, from single-energy-group (or two temperature) such as VISPHOT ( Ensman fc Bur 



rows 1992|) or TITAN ( |Gehmyer fc Mihalas 1994]) , to multigroup radiative transfer, such 



as STELLA ( Blinnikov fc Bartunov 1993 ). Other authors built codes devoted more specifi- 



cally to the radiation transport, but at the expense of detailed hydrodynamics. An example 
is the code EDDINGTON (Eastman fc Pinto 1993|) in which free expansion is assumed. 
The hydrodynamics part of most of the existing codes is a one dimensional implicit finite 
difference scheme, including artificial viscosity terms for problems requiring an accurate 
treatment of shock waves and discontinuities. During the last decade, a new subclass of 
numerical methods, the so-called Godunov-type methods, has been gradually substituting 
the schemes based on numerical viscosity due to their easier extension to multidimensional 

* We will use the term radiation for both photons and neutrinos 
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Hyperbolic character of the angular moment equations 3 
cases and their greater capabilities in the treatment of strong shocks. The lack of a radiation 

transport method fully compatible with hydrodynamical Godunov-type schemes is one of 
the motivations of these paper. 

Turning back to the radiation transport part, another important issue related to the 
angular moment equations is the closure relation. By taking angular moments of the dis- 
tribution function the complexity of BE is highly reduced, but the information about the 
angular dependence of the radiation field is partially lost. Each n th angular moment equa- 
tion of the BE contains angular moments of higher order, thus any truncated hierarchy of 
moment equations contains more unknowns than equations and must be supplied with ad- 
ditional equations or closure relations flCernohorsky fc Bludman 1994j ; proth et al. 19961 ). 



Theoretically, if the correct closure for a given problem is known, the solution obtained for 
the first moments of the distribution function by solving the moment equations should be 
the same as the solution obtained by solving the BE. 

Although the moment equations are much simpler to solve than the BE, in many situ- 
ations an additional simplification is made by neglecting some other terms, leading to the 
diffusion approximation (DA). In this approximation, however, the resulting energy flux can 
be higher than the limit predicted by causality, specially in the regions where the radia- 
tion mean free path becomes comparable to the characteristic length. Different extensions of 
DA, like flux-limited diffusion ( Powers fc Wilson 1982|) or artificial opacities ( Pgani fc Janka 



1992| ) have being used to overcome this problem. These extensions only partially solve the 
breakdown of the DA in the semi-transparent and transparent layers, and all of them are 
based on the same idea: to reintroduce some terms which had been dropped in the original 
assumptions of the DA. The term which is usually kept out of the equations is the time 
derivative of the fluxes. By neglecting this term, the character of the system of equations is 
changed to parabolic, and causality is therefore violated (disturbances propagate at infinite 
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velocity). In order to develop a method which is valid in all regimes (from diffusion to free 
streaming) while preserving causality, one must solve the full set of equations keeping its 
hyperbolic character. 

In this paper, we address the problem of establishing a well-defined hyperbolic system 
of equations for the first angular moments of the BE in the non-stationary spherically sym- 
metric case. We will see that, besides the important question of its hyperbolicity, some 
constraints on the mathematical properties of the closure relations can be derived. They 
might help to disregard among some of the closures previously proposed. Moreover, the be- 
haviour of the characteristic fields of the hyperbolic system, which gives information on the 
velocity of perturbations, can also be used to establish additional constraints to the closures. 

A direct consequence of considering the moment equations in hyperbolic form is that 
it permits the application of powerful numerical techniques that have been developed in 
recent years for hyperbolic systems of conservation laws (the equations of classical and 
relativistic fluid dynamics for perfect fluids, for example). Among the different numerical 
techniques, the so-called high resolution shock capturing (HRSC) methods have a number of 
interesting features such that stability, being conservative, convergence to physical solutions 
and high accuracy in regions where the solution is smooth. HRSC methods are based on the 
resolution of local Riemann problems (an initial value problem with discontinuous data) at 
the interfaces of numerical shells, ensuring a consistent treatment of discontinuities and steep 
gradients ( podunov 1959| ). Their special relativistic extension ( [Marti, Ibanez fc Miralles 



1991| ) has shown its potential in simulations of heavy ion collisions and extra-galactic jets, 



and different attempts to extend the method to General Relativity have been done ( Banyuls 



et ai. 199^ ; |Pons et al. 1998| ). We refer the interested reader to the recent reviews by Marti 



& Miiller ( |1999| ) and Ibanez & Marti ( |1999|) for a detailed description of the current status 



of HRSC techniques in numerical relativistic hydrodynamics. 
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Hyperbolic character of the angular moment equations 5 
Although HRSC are specially designed for hyperbolic systems without source terms 

which corresponds to the transparent regime, we will show in numerical experiments that, 
by appropriate treatment of the sources, accurate solutions can be obtained also in the 
diffusion regime, where the source terms are dominant. 

The structure of the paper is the following: In §2 we summarise the deduction of the angu- 
lar moment equations of the BE in a spherically symmetric case, and for a static background. 
We also describe the general form of the collision terms when emission-absorption and iso- 
energetic scattering processes are included. In §3 we discuss the most common techniques 
used to solve the transport equations, such as diffusion or flux-limited diffusion, remarking 
their main features and limitations. In §4 the hyperbolic character of the equations and its 
implications for the closure relations are discussed. In §5, the numerical techniques employed 
to solve the transport equations as a hyperbolic system of conservation laws are discussed. 
A number of numerical experiments solving the transport equations in several test problems 
are displayed in §6, and the feasibility of the hyperbolic treatment in all kinds of regimes is 
demonstrated. Finally, main conclusions and advantages of our proposal are summarised in 
§7. 

2 TRANSPORT EQUATIONS IN SPHERICAL COORDINATES 

We shall start our discussion from the radiative transfer equations in a static medium, 
deferring to a future work the inclusion of fluid velocity terms in the equations. Although 
the restriction to zero velocity seems to be specific, there are some astrophysical scenarios 
where the assumption that matter is at rest is a reasonable approach. Moreover, the inclusion 
of velocity terms, in some cases of interest does not change the essentials of our conclusions. 
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In a static background, the Boltzmann Equation for massless particles in spherical coor- 
dinates can be written as follows^] 



dl + dl | (l-fj, 2 ) dl 
dt ^ dr r d\i 



'dT 

dr 



(1) 

coll 



where X = X(t,r,u, /i) is the invariant distribution function, /z is the cosine of the angle of 
the particle momentum with respect to the radial direction, u the energy of the particle and 
the right hand side term is the invariant source or collision term coming from the interaction 
between the radiation field and matter. 

As stated in the introduction, a common method of solving the BE is the method of 
moments ( |Thorne 198 1|) , which involves taking angular moments of the equation by applying 



the operator 



, u l dud<& 1 = 0,1,2,... (2) 

47T J-l JO 

Denoting by E, F and P the angular moments of the specific intensity guj 3 X, g being 
the statistical weight (g = 2 for photons and g = 1 for neutrinos) 
UJ 3 [+ 1 

E = E(t, r, u) = g— J dfx X, 
on 3 [+ 1 

F = F(t, r, u) = g— J d/i fj, J, 

P = P(t,r,u)=g^ dfi fi 2 J, (3) 
the equations corresponding to the first two moments (0 and 1) of the BE are: 



IF 

d t E + d r F + — = s° (4) 
r 

d t F + d r P + ^—^ = s 1 (5) 

r 

where the moments of the collision term are defined as 

t We work in units where c= h = 1. 
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Hyperbolic character of the angular moment equations 7 

uj 2 [+ 1 (dX\ , 

Notice that, in the absence of velocity fields or strong gravitational fields, the differ- 
ent energy groups in a multifrequancy scheme are only coupled through the source terms, 
rendering the mathematical character of the left hand side of equations P])-(PD formally 
identical to the energy-averaged problem. 

2.1 Collision terms 

The most general form for the collision term including emission, absorption and isoenergetic 
scattering is the following: 

| = UJ (j - K a l + K s [l]) (8) 
coll 

where j is the emissivity, n a the absorption opacity including final states blocking (for 
fermions) or stimulated emission (for bosons) and k s the scattering opacity related to the 
scattering reaction rate (R s ) through 

K a (u)=[ — ) / dp! \ dip R s (uj, cos 6) {2{uj, //) - T(u, //)) (9) 
\2nJ J-i Jo 

where is the angle between the in-going and outgoing particle and ip is the azimuthal 
angle. The opacities and the emissivity are expressed in units of inverse length. 

For isotropic scattering the source terms appearing in equations (f|) and (B|) can be 
written as 

s° = K a {E^-E) (10) 

s 1 = -kF (11) 

where k = n a + k s and E eq is the value of E in equilibrium with matter. In the above, we 
have assumed matter in local thermodynamic equilibrium (LTE), and only the radiation 
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field (photons or neutrinos) is allowed to deviate from equilibrium. The conclusions of our 
study might be affected by large velocity fields or for frequencies close to line discontinuities, 
which are present in some astrophysical scenarios, and deserve a more careful study (|Mihalas 
et al. 197^ ; |Kunasz 1983| ; [Mihalas fc Kunasz 1986| ). 



3 FLUX-LIMITED DIFFUSION AND ARTIFICIAL OPACITIES 

One of the approaches most widely used to solve the transport equations numerically is the 
diffusion approximation, in which the invariant distribution function is assumed to be nearly 
isotropic in the comoving frame. Thus an expansion in terms of Legendre polynomials to 
the order 0(/j, 2 ) is enough to maintain the main features of the radiation field 

= J + 3 2i n (12) 

where X x <I . Consistent with this assumption, the time derivative of the flux is neglected 
in equation (|5]), and together with equation (|7]) gives the following relation for the flux in 
terms of the energy gradient 

m = -s-L^ (is, 

6K{UJ) or 

As stated before, the DA breaks down when the mean free path is large compared to 
the typical scale of the problem, and the fluxes calculated with the former formula may 
give non-causal values, in the sense that the flux can be greater than the energy density. 
This pathology comes from the fact that we have neglected some terms in the momentum 
equation in order to obtain a simple formula for the fluxes. 

Flux limiters have been introduced ad hoc in the literature to avoid this non-causal 
behaviour; see, e.g., Minerbo (1978), Levermore & Pomraining ( J1981 ) or Cernohorsky & 



Bludman (|1994j ). In order to illustrate where the problems stem from, we begin by deriving 
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Hyperbolic character of the angular moment equations 9 
the flux-limited diffusion equations. Let us define the flux factor, f, and the Eddington factor, 

p: 

/=f (M) 
P=f (15) 

By subtracting equation @ multiplied by / from (|5|), and after some algebra (now keeping 
all terms), one obtains: 

F = - } P ~ f2 \ d r E (16) 

[K ph + Kj) 

where the quantities, including the flux and Eddington factors, are now energy dependent. 
The explicit expression for the different opacities are 
E eq 

Kph = K s + K a — (17) 
fKj = d t f + d rP - fd r f + - 1 ~ (18) 

r 

The physical opacity, n p h, or effective albedo only depends on the interaction of radiation 
with matter and includes an out-of-equilibrium correction to absorption opacity. The second 
term in the denominator, kj, is the multigroup extension of the so-called artificial opacity 
( Pgani fc Janka 1992| ) and takes into account geometrical corrections and deviations from 



the diffusion limit (it vanishes in the limit /< 1). 

The term in the numerator of equation fll6|), (p — f 2 ), plays the role of a flux-limiter. It 
has the value | in the diffusion limit and goes to zero as the free streaming limit (/ — > 1 , 
p — > 1) is reached. To close the set of equations, a closure relation for p is needed, and it is 
directly related to flux limiter (p — f 2 ). We will re-address the closure relation problem in 
next section. 

With this description of the transport equations, the simple structure of the diffusion 
equation and its parabolic character is maintained when the first term in kj is neglected 
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10 J. A. Pons, J. M— Ibdnez and J. A. Miralles 

(d t f = 0). The geometric and velocity field corrections can also be included through artificial 
opacities as corrections to the physical opacity Although the flux-limiter helps to avoid non- 
causal behaviour in the sense F > E, causality is also violated when the time derivative 
of / is neglected, because the character of the equations change to parabolic, which results 
in an infinite speed of propagation of the information. By keeping this term, the diffusion 
scheme is, in some applications, numerically unstable and further numerical tools are needed 
to handle these instabilities. This is the subject of the next section. 



4 HYPERBOLIC FORMULATION OF TRANSPORT EQUATIONS 

The system of equations (Q) and ([5]) can be written in a compact form as follows 

dU_ 1 d (r 2 f) _ 
dt r 2 dr 

where the vectors U and T are 

U=(E,F) (20) 
f=(F,pE) (21) 
and the sources S, according to equations C p70| - p7Tf ) , are given by 



S = (K a {E e « - E), -kF + ^j^-) ■ (22) 
Notice that the sources are free of partial derivative operators. 



The above system of equations (19 ) is said to be a hyperbolic system of conservation 
laws (with source terms) when the Jacobian matrix associated to the fluxes ,V^JF, has real 
eigenvalues and there exists a complete set of corresponding eigenvectors. To analyse the 
Jacobian matrix, we will assume that the Eddington factor is only a function of the flux 
factor p = p(f). This assumption is commonly used by a number of authors who have 
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Hyperbolic character of the angular moment equations 11 
proposed different closure relations or flux-limiters in the context of flux-limited diffusion. 



In this case, the Jacobian matrix is 
/ 1 



V p- fp' p' , 

dp 



(23) 



where p' . 

df 

The eigenvalues of the above matrix are 



p' ± yV 2 + 4(p - fp>) 

A ± = 2 { ' 

and the right eigenvectors are 

f + = (l,A + ); rL=(l,A_) (25) 

The mathematical character of the system of equations ( |19|) depends on the value of the 
discriminant in equation ([24]), which can be written as follows 

A = (p'-2f) 2 + 4(p-f 2 ). (26) 

Since p and / are normalised moments of a non-negative weight function, on the unit sphere, 
there is an extra relation between / and p, the Schwarz inequality 

f<P<l , (27) 

that ensures that for any closure satisfying ( pT| ) the eigenvalues are always real and hence 
the hyperbolic character is guaranteed. In addition, causality imposes that the velocity of 
propagation of perturbations must be smaller than the speed of light, which gives another 
constraint on the eigenvalues 

|A±| < 1 , (28) 

that can be expressed in terms of /, p and p' as follows 

<!/<?-*. (29) 
1 + / ~ P ~ l-f { ' 
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12 J. A. Pons, J. M— Ibdnez and J. A. Miralles 

This restriction, together with equation ( p7[ ) and the condition |/| < 1, defines the region 
in the [/, p, p'] space allowed for causal closures. 

In many astrophysical problems, radiation diffuses out from a central opaque region to an 
outer optically thiner region. In such a scenario, the radiation angular distribution is nearly 
isotropic in the central region (|/| <C 1) and forward peaked in the transparent region, far 
from the opaque centre (/ — > 1). Thus, the closure in that problem must satisfy 

p(/ = 0) = i and p(/ = l) = l. (30) 

The first limit comes from the diffusion approximation while the second limit comes from 
the Schwarz inequality (p7|), and both are satisfied by all the closures found in the literature 
(see Table I). Another general property is that, in the diffusion regime the Eddington factor 
has the constant value of 1/3, thus another restriction on the derivative can be imposed 

P'(f = 0) = , (31) 

that leads to the characteristic speeds A± = ±i. In the limit / = 1, equations fl27|) and ( 29|) 
gives the bounds < p 1 < 2 and the eigenvalues in this limit are A + = 1 and A_ = p 1 — 1. 
According to the theory of hyperbolic systems, this result can be interpreted as a first wave 
propagating outwards at the speed of light and a second wave with a characteristic speed 
that depends on the value of p'. Although A_ depends on the form of the closure, for / = 1 
the amplitude of this second wave is zero and the value of p' does not affect the evolution. 
However, for / close but not equal to one, the choice of p' can be important, as we will show 
in section §5.2. Similar conclusions can be obtained for the case / = —1. 

In table 1 we summarise the main properties of some of the most widely used closures: 
Minerbo (1978), LP ( |Levermore &; Bomraining 1981 ), the linear formula by Auer ( |1984 ). 



Kershaw's parabolic approach (Kershaw 1976|) , the flux-limiters used by Bruenn ( 1985[) 



and Bowers & Wilson (|1982|) , Janka's fit to Montecarlo simulations of neutrino transport 
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Hyperbolic character of the angular moment equations 13 



in Supernovae flJanka 199 1| ) and that corresponding to an isotropically emitting sphere 



flUoopcrstcin, van den Horn fc Baron 1986|) , named vacuum. 



There are several comments about Table I that we would like to address. First, one closure 
relation (BW) leads to non-causal values of one eigenvalue in the diffusion regime. Secondly, 
four closures do not satisfy the condition p'(f = 0) = (Auer, Bruenn, BW, Vacuum). Auer's 
formula failure is due to an excessive simplification (it is a linear approach) while 'Vacuum' 
has been obtained from the assumption of an isotropically radiating sphere, which restricts 
their validity to / > 0.5. The closure relations proposed by Kershaw, Janka, Minerbo and 
LP have the common features of preserving causality leading to reasonable values of the 
characteristic velocities (A± = ±—7^) in the diffusion regime. Thirdly, it is worthy to point 
out that all the closures except Auer's give p' > 1 when approaching the radially-streaming 
limit (/ = 1), which makes both eigenvalues be positive. It ensures that no information can 
be propagated inwards for sufficiently high values of /. It is clear that these closures cannot 
be used to treat problems in which a small perturbation is propagating inwards from the 
transparent region, as we will show in next section. 

In a more general situation, the Eddington factor depends on both, the flux factor and the 
occupation density corresponding to the zeroth angular moment of the distribution function, 
e = E/u 3 . The Jacobian matrix associated to the fluxes is 







1 \ 



dU 




and now the eigenvalues are 



dp 

die 




(32) 



where the discriminant is 
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(33) 

In contrast to (PUD, the sign of the discriminant might be negative for negative values 
of |^ , and the system would become elliptic instead of hyperbolic. As stated before, hy- 
perbolic equations are related to finite speeds of propagation of the perturbations, while 
elliptic equations are related with boundary problems and do not allow for solutions with 
propagating signals. This fact indicates that, in order to preserve hyperbolicity, only closures 
which give positive values of the discriminant should be considered. 

A closure of this kind has been proposed by Cernohorsky & Bludman ( |1994|) using ar- 



guments of maximum entropy. Minerbo's, LP's and the Vacuum closures can be obtained 
as limiting cases of CB closure in the low occupation limit (classical statistics), large ocup- 
pation limit for Bose-Einstein statistics and maximal forward packing limit for Fermi-Dirac 
statistics, respectively. We have checked that CB closure gives real and causal eigenvalues 
for both, Bose-Einstein and Fermi-Dirac statistics, in any regime. 

5 THE NUMERICAL APPROACH. 

The numerical code employed for the calculations presented in §6 is a finite finite difference 
scheme in which the left (L) and right (R) states for the Riemann problems set up at each 
interface have been obtained with a monotonic, piecewise linear reconstruction procedure 
( |van Leer 1979| ). The numerical fluxes are evaluated following the idea of Harten, Lax & van 



Leer (|1983), that is based on approximate solution of the original Riemann problem with a 



single intermediate state. The resulting numerical flux is given by 
where ip + = max(0, A+, A+) and ip~ = min(0, X R , \ L ). 
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Hyperbolic character of the angular moment equations 15 
The time-advance algorithm is a TVD-preserving, third order Runge-Kutta that takes 

into account the influence of the source terms in intermediate steps. The source term in this 
problem deserves a special mention. It acts as a relaxation term which leads to a long-time 
behaviour governed by the reduced parabolic system described in section §3. Although high 
resolution methods for hyperbolic conservation laws usually fail to capture this feature (un- 
less a very fine grid is used), the use of a modified flux (|Jin fc Levermore 19961 ) and an 



implicit linearization of the source term in each Runge-Kutta step allows us to succeed in 
the treatment of such stiff relaxation term with coarse grids and time-steps longer than the 
source time-scale (« as we will show in subsection §6.1. Although we have constrained 
the numerical experiments to situations when only scattering processes are permitted, the 
inclusion of absorption-emission processes can be treated in a similar way. However it in- 
troduces energy and momentum exchange with the matter and makes necessary to solve 
the equations of radiation transport coupled with the equations of hydrodynamics or the 
equations of hydrostatic equilibrium with number and heat transfer. Such a problem is out 
of the scope of this paper, in which we focussed on the radiation part, and will be addressed 
in forthcoming works. 

6 NUMERICAL EXPERIMENTS AND DISCUSSION. 

To check the applicability of HRSC techniques to radiation transport problems, we have 
performed a number of numerical tests considering several simplified models in which the 
analytical solution is known. In the remainder of this section we describe the different tests 
and display the results. Our intention is to prove that our method can solve radiation 
transport in all kinds of regimes, from diffusion (TEST 1) to radially streaming (TEST 
2), going through the semi-transparent region (TEST 5) and in situations where strong 
gradients are developed (TEST 4). 
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Ideally, in a multi-group scheme we have as many systems of two equations as energy 

bins in which the radiation density is splitted. In a radiation hydrodynamics problem all 

those pairs of equations are coupled to the matter equations and must be solved together. 

Since our objective is to study the applicability of HRSC techniques for the radiation part, 

in what follows we assume to have only one energy group or, equivalently, we consider the 

energy integrated equations. Notice, though, that the conclusions are not affected since, 

when velocity fields are not taken into account, the equations for different energies are only 

coupled through the source terms. A different problem will arise when velocity fields are 

important or gravitational and Doppler red-shift are considered. In those situations either 

the anisotropy of the source term (in the observer frame) or the advection and aberration 

terms (in the comoving frame) cause the energy redistribution along the photon (or neutrino) 

path and a more careful treatment needs to be done. 

6.1 TEST 1: Diffusion limit 

In the diffusion limit, the closure is obviously p — 1/3 and the time derivative of the flux 
can be neglected, leading to equation (|l|). Setting the absorption opacity to zero and the 
scattering opacity to a constant value, the analytical solution of the system formed by (£|) 
and (|13D , corresponding to an initial value problem consisting in a Dirac delta located in 
r = at t = 0, is the following 



We imposed F(r = 0) = as inner boundary conditions and outflow boundary conditions 
in the outermost shell. 




(35) 




(36) 
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Hyperbolic character of the angular moment equations 17 



We have performed two simulations. First, we set a value of k = 100 on a equally spaced 
grid of 100 shells in the domain [0,1]. It corresponds to a Peclet number (Pe = kAx) of 
order unity. The cell Peclet number is a measure of the spatial resolution relative to the 
relaxation scale and a grid is coarse when Pe ^> 1. This first simulation starts at t — 1 to 
avoid the numerical problems produced by an infinite value in t = 0. Results are shown in 
Figure [l|. The energy density (left panel) and the flux (right panel) are displayed for different 
times in dimension-less units, t = 1, 2, 3, 5 from top to bottom. The solid line corresponds 
to the analytical solution and the crosses correspond to the numerical solution. Notice that, 
although we are finding numerical solutions of system (|19"D, the analytical solution of the 
diffusion limit is well reproduced by the numerical code, since the opacity is large enough 
to be close to the diffusion limit. The second simulation has been performed in a grid of 100 
shells in the domain [0, 0.5], and using a value of k = 10 5 , which corresponds to Pe = 5 x 10 4 . 
The initial time in this case is 200 and the Courant time-step is of the same order (10 -2 ) 
as in the previous simulation, three orders of magnitude larger than the relaxation scale. 
Results are displayed in Fig. ^|. The agreement between the analytical and the numerical 
solutions clearly proofs the capability of the code to handle situations with very stiff source 
terms. 

6.2 TEST 2: Radially Streaming limit 

In the radially streaming limit (|/| = 1, p = 1), and assuming there is no interaction between 
radiation and matter (/t a = 0, k — 0) equations ( |T9"D are reduced to 




uj) 1 d(r 







(37) 



Or 
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The solution to this equation is any linear combination of ingoing and outgoing spherical 
waves, i.e., any linear combination of functions of the form E(t, r) = sg(r ±t). For the first 
numerical experiment we choose a Gaussian spherical wave propagating outwards, 

E(r,t) = -^exp{-(r-t) 2 } (38) 

F(r,t) = E(r,t). (39) 

Our numerical grid consists on 200 equally spaced shells going from r\ = 0.2 to r2 = 10.2 
and the boundary conditions consist in imposing the analytical solution in both boundaries. 
We show the results of this test in Figure |3|. In the figure we displayed the analytical solution 
(solid line) and the numerical (crosses) for the energy density. In order to analyse the effect 
of the particular choice for the closure, we repeated the experiment using different values for 
p'(f = 1), which leads to different speeds of propagation of the waves, and no difference was 
observed. The final reason is that all the closures displayed in Table 1 give A+(/ = 1) = 1. 
This is a general property of any closure, since p — f 2 is a positive defined function that goes 
to zero as / — > 1 and, therefore, its derivative has to be negative in a vicinity of / = 1 and 
p' < 2f. Thus taking the limit / — ► 1 in (|24[) , one obtains 



p' + J(2-p') 2 

A+ = % = 1 (40) 

To understand why the value of A_ is not important in this problem, we study the 
diagonalized system obtained by adding and subtracting the two equations in ( |T9"D after 
taking p — 1, 

at r 2 or 

at r 2 or 

In this equations it is clearly seen that the quantities propagated with velocities A± (Riemann 
invariants) are E ± F, respectively. Since in our particular problem we took F = E, the 
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quantity propagated by A_ is identically null and the value of the speed of propagation does 

not affect the solution. 

A different situation is found when we add to the initial conditions a second wave propa- 
gating in the opposite direction. Figure [| shows the results from such experiment, in which a 
second wave with smaller amplitude and opposite direction is located initially at r = 7. The 
closure used was p = 1, p' = 0, which gives A + = 1, A_ = — 1. In the figure, the analytical 
solution is denoted by the solid lines and the numerical solution by crosses and each panel 
corresponds to a different time. It can be seen how the code is able to solve the interaction 
of the two waves until they completely cross over. The difference between the analytical and 
numerical solution near the center for the last panel are due to the boundary conditions 
(initially, it was an ingoing wave). 

To illustrate the effect of a wrong speed of propagation, the experiment is repeated with 
Kershaw's closure, that gives a correct speed of propagation when / = ±1 but a wrong 
value otherwise. Results are shown in Figure f| for the same time steps as in the previous 
case. At the beginning both waves seem to propagate correctly but, as they cross each 
other and the flux factor departs from unity, the numerical result differs from the analytical 
solution. A similar behaviour is obtained with other closures. This simple example manifests 
the importance of the choice for the closure, since it fixes the speeds of propagation of the 
waves. 

6.3 TEST 3: Sphere radiating in the vacuum. 

In the two previous tests we studied the limiting cases, diffusive and free-streaming regimes. 
Another academic problem in the limit of vanishing opacities consists in a sphere of radius 
R radiating isotropically into vacuum (Figure |6|) . Given the explicit form of the distribution 
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function at the surface of the sphere, Z (£) = X{t = 0, r 

easily obtained 

(lo(t-s) n>x{r,t) 
I(t, r,fi) = l 

[ fji<x(r,t) 

where 



R), the general solution can be 



(43) 



s = r 



x(r, t) 



i V 1 - (w 

t 2 + r 2 - # 2 
2tr 

1 



t > Vr 2 - R 2 



r - R<t< Vr 2 - # 2 
t<r-R 



(44) 



(45) 



where /x = cos a. 

We take a linear time dependence Zo(i) = Ct, being C an arbitrary constant, so that the 
angular integrals can be done analytically The angular moments of the distribution function 
are, then, given by 
C 



E(t,r) 
F(t,r) 
P(t,r) 



4 

£ r 

12 

£ 

48 



2t(l — x) — rfl — a; 2 ) + rVl — x 2 + rx 2 log [ X 

\l + ^/T^x 2 



3t(l - a; 2 ) - 2r(l - x 3 ) + 2r(l - x 2 f /2 



- x 3 ) - 6r(l - a; 4 ) + 3r(2 - x 2 ) v / T 3 x 2 + 3rx 4 log 



x 



(46) 
(47) 
.(48) 



.1 + Vl-x 2 , 

We set up a numerical grid of 200 points, uniformly distributed, from r = 1 to r = 200. 
The initial conditions are F = and E = a/r 2 with a being a small number, different from 
zero for numerical purposes. We have taken a = 10~ 4 . The radiation field is introduced 
through the inner boundary conditions, E(t,r = R) = 2t and f(t,r = R) = 0.5. Outflow is 
assumed in the outer boundary. 

In Figure [7] the evolution of the radiation energy density as a function of the radial 
coordinate is shown for three different evolution times. The test has been repeat for three 
different closure relations from those described in table 1: Vacuum (dashed lines), Kershaw 
(dotted lines) and Minerbo (dashed-dotted lines). As in the previous figures the analytical 
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solution is represented by the solid line. Although the appropriate closure for this problem 
(vacuum) gives a better approximation to the real solution, it is remarkable that other 
closures also give similar results. The differences between different closures stem from the 
slightly different propagation speeds of the simple waves in the regime of interest. Unlike in 
the two-wave problem described in previous section, for this specific problem any reasonable 
closure might be used without major deviations from the real solution. 

6.4 TEST 4: Radiative cooling and heating 

With this test we intend to simulate a more realistic scenario in which, initially, radiation 
produced in the center of a star diffuses out with constant luminosity. We begin by setting up 
a sphere which is in radiative equilibrium at a given luminosity L. For p—\, the stationary 
solution is given by 



where R is the radius of the outer boundary, at which / = 1. Starting with the initial model 
corresponding to L — 4ir, we let it evolve in three different situations: 1) we take L = as 
boundary conditions at the inner boundary (located at r = 1), and the radiation field should 
diffuse out continuously; 2) we change the inner boundary condition to L — 40ir, ten times 
larger than the initial model. The behaviour is now the opposite, the energy density must 
increase until the new stationary solution for the new luminosity is reached; 3) we again take 
L = at the inner boundary but we switch on a central point-source in the middle of the 
star (r = 6). The radiation field must evolve towards the corresponding stationary solution, 
which is different in the two parts of the star defined by the source. Inside, the luminosity 
is zero and the energy density is constant. In the region external to the energy source, the 
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(49) 



F(r) 



L 



(50) 



4nr 2 
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stationary solution is denned by L = s, where s is the energy per unit time injected by the 
source. The results for k = 100 are shown in Figures |8|-|T0|, for cases 1 to 3, respectively. 
All three simulations have been done on the same grid of 200 uniform zones and outflow 
boundary conditions in the outer edge have been imposed. 

Left panels display the energy density as a function of radius. The thick solid line in the 
figures corresponds to the initial model and the thick dotted line in Figure |9| represent the 
stationary solution with L = 40ir. In Figure ||, we can see how the heat wave lasts about 
1000 time-units (the diffusion time Rk) to arrive to the surface and a new stationary solution 
is successfully obtained after m 15 times the diffusion time scale. 

6.5 TEST 5: The semi-transparent regime. 

Testing a method in the semi-transparent regime becomes a tough issue due to a lack of 
analytical solutions to compare with. To overcome this problem, we have performed Monte 
Carlo simulations, which is the closest solution to an exact one, and results in an excellent 
test of the performance of the method in the semi-transparent regime. 

We have used a Monte Carlo code that simulates the transfer of a radiation field in a 
spherically symmetric geometry. We consider a region with inner boundary at r = R in and 
outer boundary at r = R ou t, divided in 200 spherical shells. Within each shell a constant 
value of the scattering opacity (elastic and isotropic) is assumed. We do not consider emission 
or absorption processes. Our Monte Carlo procedure starts by injecting particles at the 
inner boundary outwardly-directed. Then, we compute the trajectory of each particle as it 
scatters off matter, according to the standard random walk laws ( [Lucy 1999| ). The trajectory 



ends when the particle escapes by crossing the outer boundary or it re-enters the inner 
boundary. Once the path of a particle is computed we proceed to calculate the contribution 
of this particle to the moments of the radiation field. The contribution of each particle 
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to the energy density at a given numerical shell is calculated by adding the path length 

between two successive scattering events, within the shell, divided by the volume of the 

shell. Consequently, the total energy density is obtained just by adding the contribution of 

each particle multiplied by an arbitrary factor Analogously, to obtain the contribution of 

each particle to the flux, the path length is weighted using the cosine of the angle with respect 

to the radial direction. Finally, pressure is calculated by means of the same procedure, but 

using the square of the cosine as weight of the path length. 

In our numerical experiment we have taken: R in = 10 and R out = 100. The opacity law is 
an exponentially decreasing function in radius, taking the value k = 2 at the inner boundary 
and k = 2 x 10~ 5 at the outer boundary. The corresponding values of the optical depth are 
such that an important part of the region is semi-transparent. We have run the Monte Carlo 
code using 500 000 particles in order to obtain the stationary solution. 

The set up for our transport code is the following: i) the above opacity law, ii) the inner 
boundary conditions given by the stationary solution (obtained from Monte Carlo), and iii) 
the closure relation. We have performed two different calculations using Kershaw's closure 
relation and a cubic polynomial fit to Monte Carlo results, explicitly 

p = 0.3228 + 0.1902/ - 0.0476/ 2 + 0.5131/ 3 (51) 

Starting with arbitrary initial profiles, we let the radiation field evolve until the stationary 
solution is reached and results can be compared to those obtained from the Monte Carlo 
simulations. 

In Figure [Li] we show the stationary solutions calculated from Monte Carlo simulations 
(solid line) and from our code with the two different closures, Kershaw's (dotted line) and 
the above fit (dashed line). From the results displayed in the figure we can conclude that the 
evolutionary code describes the semi-transparent regime accurately. The differences in the 
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values of the energy density (left panel), normalized to its inner boundary value, amounts 
less than a few percentage in any case. The slight discrepancies in the flux factor are mainly 
due to the closure relation employed, being less than 6% for Kershaw's closure and less than 
3% for the fit. The closure, therefore, sets the maximum accuracy that can be obtained, 
although qualitative results at a level of 10% uncertainty can easily be obtained with any 
reasonable closure in the semi-transparent regime. 

6.6 A toy model. 

Our last experiment intends to illustrate how the numerical method we proposed can succeed 

in situations when standard flux-limited diffusion techniques fail. We set up a sphere of radius 

r = 5 with a radiation energy density of E = 1, embedded in a background where the opacity 

varies as follows: 

k r < 1 

K /r 2 1 < r < 10 

/.' < 10 < r < 15 (52) 
k 15 < r < 17 
r > 17 

where we have taken k = 10 3 . Outside the sphere we take E — (in practice, E = 10~ 6 
for numerical reasons), and initially the flux is zero everywhere. It tries to mimic a situation 
in which radiation scatters out in a central star, where the diffusion approximation remains 
valid, through an opacity- varying atmosphere up to the surface (r = 10), where the opacity 
suddenly vanishes. Then, it crosses a transparent region to find a very opaque layer, where 
diffusion operates again. This sort of problem is obviously a tough challenge for diffusion- 
based codes, since they cannot work properly in the transparent region. 



In Figure we show results from the evolution of this problem on a grid of 200 uniform 
shells and using Kershaw's closure. In the left panel we show the energy density for different 
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times. It can be seen how a gradient is developed in the diffusive regions (r < 10, 15 < r < 17) 

while radiative equilibrium (dE/dr = 0) in the transparent region between the two opaque 

layers is rapidly obtained. Notice also that only for the last snapshot (t = 1000, dashed 

curve) there is a flux of energy in the outer transparent region, because it takes some time 

to the radiation to diffuse through the opaque layer and develop a gradient that keeps the 

outgoing flow. 

7 CONCLUSIONS 

We have studied the mathematical character of the angular moment equations of radiative 
transfer, considering the implications of using different closure relations. The role played 
by the closure and their derivatives is more apparent in the hyperbolic formulation. After 
this analysis, we conclude that the hyperbolic character is assured for a given variety of clo- 
sures widely used in the literature: those given by p = p(f). For general closures p = p(e, /), 
hyperbolicity can not be guaranteed, although the closure proposed by Cernohorsky & Blud- 
man fljggg) , based on maximum entropy arguments, leads to real eigenvalues. Additional 



constraints on the closures are imposed on the base of the behaviour of the eigenvalues of 
the Jacobian matrix since they give the velocity of propagation of perturbations. Causality 
limitation is written in terms of the closure relation helping us to select, from this point 
of view, among the different closures reviewed. It turns out that only one of the closures 
analysed does not satisfy causality, allowing for velocities of the waves higher than the speed 
of light. For the other closures, although derived without taking into account the mathe- 
matical issues addressed in this paper, we found that they are consistent with hyperbolicity 
and causality. 

Writing the moment equations of the distribution function as a hyperbolic system of 
conservation laws allows one to apply numerical techniques specifically designed for such 
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systems, the so-called HRSC methods. We have applied HRSC methods to solve the trans- 
port equations in a static background showing, for a number of test problems, the feasibility 
of the method. From the results, the main conclusion is that the numerical method can be 
used in any regime, from optically thick to transparent regions, obtaining numerical solu- 
tions with high accuracy. It turns out that the closure relation plays an important role, more 
evident than in traditional methods for radiation transfer, that can be useful to choose the 
more appropriate closure for a given problem. 

When the radiation transport equations are written as a system of conservation laws, 
the coupling with hydrodynamics is straightforward ( [Lowrie, Morel fc Hittinger 1999|) . This 
feature make us to be confident in the possibility of applying the method to problems 
involving radiative flows. Moreover, we are optimistic in the applicability to multidimensional 
problems, since the extension of the method is straightforward. 
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Table 1. Closure and characteristic structure. 
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Figure 1. Moderate diffusion limit: Energy density (left panel) and flux (right panel) as a function of the radial coordinate 
for different times, (t = 1, 2, 3, 5 from top to bottom). The opacity has been set to k = fOO, corresponding to a Peclet number 
of Pe = 1. 
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Figure 2. Extreme diffusion limit: Energy density (left panel) and flux (right panel) as a function of the radial coordinate for 
different times, (t = 200,240,300,400 from top to bottom). The opacity has been set to k = 10 5 , corresponding to a Pcclet 
number of Pe = 5 X 10 4 . 
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Figure 3. Free streaming limit: A Gaussian spherical wave propagating outwards. The analytical solution is represented by 
the solid line and the numerical solution by crosses. Snapshots for four different times (t = 2.5, 3.5, 5.0, 7.0) are plotted. 
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Figure 4. Free streaming limit: Two Gaussian spherical waves crossing over. The solid line corresponds to the analytical 
solution and the crosses to the numerical solution. Each panel displays a snapshot for a different time (t = 2.5,3.5,4.5,6.0). 
The closure employed is p = 1, p' = 0, which gives \± = ±1 everywhere. 
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Figure 5. Same as Figure W but using Kershaw's closure. The effect of taking wrong characteristic speeds is clearly illustrated. 
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Figure 6. Schematic diagram of Test 3: Sphere radiating isotropically into vacuum. The analytical solution for the radiation 
distribution function can be obtained in terms of r and fj, = coso. 
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Figure 7. Sphere radiating in the vacuum: Energy density as a function of the radial coordinate for three different evolution 
times. The analytical solution (described in §5.3) is represented by solid lines while the numerical solutions obtained with three 
different closures are represented by dashes (vacuum), dots (Kershaw) and the dash-dotted line (Mincrbo). 



Ld 



1000.0 



100.0 



1 0.0 



.0 



0.1 



t — i — i — | — i — i — i — i — i — i — i — | — i — i — i — | — i — i — i — i — i — i — r 




_l i i i I i i i I i i i I i i i L 



1 .2 
1 .0 

0.8 

0.6 

0.4 

0.2 
0.0 



t — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — r 



t = 500 




J I I I I I I i r I 



2 4 6 8 10 12 



2 4 6 8 10 12 



r 



r 



Figure 8. Radiation initially in radiative equilibrium (L = An) diffusing out after switching off the central energy source. The 
solid thick line represents the initial model. The energy density is plotted in the left panel and the luminosity in the right panel. 
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Figure 9. Radiative heating: starting from the same initial model (solid thick line) as in figure we increase the luminosity 
at the inner boundary to L = A0n. The energy density evolves towards the new stationary solution at higher luminosity. 
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Figure 10. Same as figure g but with a point-like energy source located at r = 6. The region inner to the source reaches the 
stationary solution at L = 0, while the outer region settles in radiative equilibrium at luminosity higher than the initial model. 
The sharp discontinuity in the luminosity is well resolved. 
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Figure 11. Semi-transparent regime: Stationary solution of the radiation field in a background with a scattering opacity 
exponentially decreasing with radius. Solid lines correspond to the Monte Carlo simulation and dotted and dashed lines 
correspond to the results obtained with the evolutionary code with Kershaw's closure and a numerical fit to Monte Carlo 
results, respectively. Left panel shows the normalized (to its inner boundary value) energy density as a function of the optical 
depth. Right panel shows the flux factor. 
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Figure 12. Evolution of the radiation field in a background with strongly varying opacity. Snapshots of the energy density (left 
panel) and flux (right panel) are shown for different times t = 0.5, 100, 250, 500, 1000. The thick curve corresponds to t = 0.5 
and the dashed curve to t = 1000. See the text for details on the opacity profile. 
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